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Abstract 

Interacting fluids, endowed with bulk viscous stresses, are discussed in a unified per- 
spective with the aim of generalizing the treatment of cosmological perturbation theory 
to the case where both fluctuating decay rates and fluctuating bulk viscosity coefficients 
are simultaneously present in the relativistic plasma. A gauge-invariant treatment of 
the qualitatively new phenomena arising in this context is provided. In a complemen- 
tary approach, faithful gauge- fixed descriptions of the gravitational and hydrodynamical 
fluctuations are developed and exploited. To deepen the interplay between bulk viscous 
stresses and fluctuating decay rates, illustrative examples are proposed and discussed 
both analytically and numerically. Particular attention is paid to the coupled evolution 
of curvature and entropy fluctuations when, in the relativistic plasma, at least one of 
the interacting fluids possesses a fluctuating bulk viscosity coefficient. It is argued that 
this class of models may be usefully employed as an effective description of the decay 
of the inflaton as well as of other phenomena involving imperfect relativistic fluids. 
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1 Formulation of the problem 



It is plausible, as observations suggest, that the large-scale temperature fluctuations de- 
tected in the microwave sky may arise, via the Sachs-Wolfe effect, from primordial curvature 
perturbations amplified during the early stages of the life of the Universe. It is therefore 
mandatory for those interested in this possibility to scrutinize in detail the evolution of cur- 
vature fluctuations when the composition of the primeval plasma cannot be parametrized 
in terms of a single, perfect, relativistic fluid. 

A class of dissipative effects that can enter the energy-momentum tensor without spoil- 
ing the isotropy of the background geometry can be parametrized in terms of a bulk viscosity 
coefficient [1] (see also [2]). In fact, in the relativistic theory of imperfect fluids, the dissipa- 
tive effects can be parametrized in terms of a shear coefficient, a heat flow coefficient and a 
bulk viscosity coefficient. When written in general coordinates [1] (see also [3]), the contri- 
butions of the shear tensor and of the heat flow (unlike the bulk viscosity) do not contribute 
to a homogeneous and isotropic background geometry of Friedmann- Robertson- Walker type 

[!]■ 

Following then Ref. [1], for a single fluid, the total energy- momentum tensor can then 
be split into a perfect contribution, denoted in the following by and into an imperfect 
contribution, denoted by AT^, i.e. 

t; = r; + at;, (i.i) 

In general coordinates, and within our set of conventions, the contribution of bulk viscous 
stresses can be written, in turn, as 

AT; = e(^-nX)v^ a , (1.2) 

where £ represents the bulk viscosity coefficient. The presence of bulk viscosity can also 
be interpreted, at the level of the background, as an effective redefinition of the enthalpy 
(i.e. (p + p)) that appears both in the spatial components of Einstein equations and in 
the background form of the covariant conservation equation. While this approach may be 
appropriate at the level of the evolution of the background space-time, it could be confusing 
if applied to the dynamics of the fluctuations when £ is allowed to have a spatial dependence. 

The presence of bulk viscosity in the relativistic plasma can influence the character of the 
cosmological singularity [3, 4] (see also [6, 5]). Regular solutions [4] (but generically unstable 
[3] and probably not past-geodesically complete) can be found if the £ has some particular 
functional dependence upon the energy density. Some of these solutions have the property 
that as t — > — oo (t being the cosmic time coordinate) the scale factor approaches a quasi-de 
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Sitter stage of accelerated expansion. These deflationary solutions [7] violate, initially, the 
strong energy condition, i.e. (p + 3V) > 0, where V = p — 3 H£ is the total pressure and 
p is the pressure arising from the perfect contribution (i.e. T- = —pSf) introduced in Eq. 
(1.1). The dominant energy condition [8] could be violated when bulk viscous stresses are 
added [7]. This second occurrence implies that the effective enthalpy, i.e. (p + V), may not 
always be positive-definite in all the classes of solutions involving a homogeneous (viscous) 
background. Boundary conditions can be imposed in such a way that for t — > +oo the 
solution approaches, asymptotically, de Sitter space-time. In the latter case the solutions 
are driven to an inflationary phase [7]. 

Barrow [7] presented detailed studies (see also [9, 10]) on inflationary Universes driven 
by a gas of strings showing that, in some regimes, this system can be dicussed in terms of 
an effective bulk viscosity coefficient. The solutions quoted in the previous paragraphs have 
been generalized along different lines. In [7] the bulk viscosity coefficient £ is taken to be 
proportional to a power of the energy density, i.e. £ oc p m . The relevance of these examples 
for the cosmic no-hair theorems has also been investigated in a series of papers [9, 10]. 
Moreover it was argued that, as one traces the Universe back in time, a point is reached 
when the expansion rate is so fast that the rate of string creation due to quantum effects 
balances the dilution of the string density due to the expansion [11]. These considerations 
lead naturally to an effective bulk viscosity coefficient that may induce an approximate, 
weak, inflationary dynamics. 

The possibility of driving a robust inflationary phase by means of bulk viscous stresses 
has been indeed criticized on various grounds (see [12, 13, 14]). In [12] it has been argued 
that the bulk viscosity effects associated with a gas of weakly interacting particles cannot 
make the pressure negative excluding any form of inflation. In [13] it has been suggested 
that hydrodynamic inflation is intrinsically problematic (see also [14] for an interesting, 
related perspective) since while the inflationary evolution requires a very small expansion 
time the hydrodynamic regime requires a very small interaction time. Therefore, to treat 
consistently the problem one would need a consistent (non-linear) generalization of the 
Israel-Stewart transport equations (see the attempt of Ref. [15] ) as well as a consistent 
model of fluid behaviour under super-rapid expansion and strongly non-linear conditions. 
The possibility of realizing a sufficiently long and robust inflationary phase is, however, not 
central for the present analysis: in the present context the possible effects of bulk viscosity 
will arise after the end of a (conventional or unconventional) inflationary phase. 

The effect of bulk viscous stresses has also been studied in the case of scalar-tensor 
theories of gravity (see for instance [16]). More recently, the possible implications of bulk 
viscosity on the dynamics of the dark energy fluid have been investigated in different frame- 
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works [17] (see also [18] and references therein). In [17] the hope is to obtain naturally an 
accelerated expansion from bulk viscous effects. 

An apparently unrelated effect is the energy exchange between different fluids whose 
mixture constitutes the global relativistic plasma that drives the evolution of the geometry. 
For instance, it is conceivable that, after inflation, the oscillating inflaton field decays into 
massless particles [19]. Then this situation can be modelled by a fluid of dusty matter 
decaying into a fluid of radiation. Multifluids also arise is the discussion of isocurvature 
modes like, for instance, the CDM-radiation isocurvature mode. In all these problems, one of 
the relevant aspects is the analysis of the evolution of curvature and entropy perturbations. 

The theory of cosmological perturbations in the presence of different fluids has a long 
history. For instance, in [19] the problem of multifluids was analysed in the absence of 
bulk viscosity. In [22] multifluid systems have been treated in the framework of a covariant 
cosmological perturbation (see [23] for an introduction and [24] for related applications). 
Entropy perturbations arising in multifluid systems have also been discussed in [25, 26, 27, 
28] in the context of the Bardeen formalism [29], as well as in the context of inflationary 
models [31, 35] and quintessence models [32, 33, 34, 36]. Possible effects connected with the 
decay of the inflaton have recently been considered in [37, 38]. Notice finally that earlier 
works in the general area of bulk viscosity can be found in [39, 40]. 

It has been argued, in recent years, that interacting (multiple) fluids may play some 
relevant role in cosmological perturbation theory. If the early Universe contains more than 
a single fluid, the spatial variations of the decay rate may indeed generate curvature per- 
turbations [41] (see also [42] and [43, 44, 45, 46]). 

In the present paper a unified description of imperfect fluids is developed with the aim 
of understanding if and how the effect of spatial fluctuations in the bulk viscosity coefficient 
(combined with a non- vanishing decay rate) may induce super-horizon curvature perturba- 
tions. In this sense, we are not interested in discussing the evolution of the fluctuations in 
the case when a single viscous fluid is present. Let us phrase the problem in physical terms. 
Consider the situation where the matter content of the Universe is formed by, at least, 
two fluids. The possible dissipative effects compatible with the isotropy of the background 
geometry are: 

• energy exchange between the fluids; 

• bulk viscous effects associated with each fluid and with the mixture. 

The first type of efffect can be parametrized by a (time-dependent) decay rate that is 
allowed to fluctuate spatially. These fluctuations may then induce fluctuations on super- 
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horizon scales (see, again, [41, 42, 43, 44, 45, 46]). The second type of effects can be 
parametrized by the appropriate second viscosity coefficients of each fluid of the mixture and 
by their related fluctuations. Also in this case it is reasonable to expect that the fluctuations 
of the second viscosity coefficient can be converted into super-horizon curvature fluctuations. 
There may be, in principle, two extreme situations. In one case the presence of bulk viscous 
stresses crucially determines the nature of the background geometry. This is the approach 
studied, for instance, in [4, 7, 11, 16]. In the other situation, the bulk viscosity does not 
determine crucially the properties of the background geometry. This approach is hard to 
justify for a single (non decaying) viscous fluid, but it becomes plausible when there are at 
least two interacting fluids. Consider, for instance, the situation where the fluid describing 
the decay products has no intrinsic viscosity, while the fluid that is effectively decaying has 
some intrinsic viscosity. If, as in [7], £ depends on the energy density of the decaying fluid, 
it reasonable to expect that when the decay is completed, the bulk viscous effects will also 
decay. Depending on the specific parameters of the example, it is then conceivable that, the 
net result of the evolution will be a single inviscid fluid supplemented by an exponentially 
decaying (viscous) component with negligible effect on the dynamics of the background. 
It is interesting to understand what happens to the dynamics of the fluctuation. If the 
decaying fluid has some intrinsic viscosity, the effect on the curvature fluctuations may be 
relevant even if, when the decay is completed, the viscous contribution goes exponentially 
to zero at a rate that depends upon the relative interaction of the fluids of the mixture. 

From a more formal point of view, one of the problems is a full gauge-invariant treatment 
of imperfect cosmological perturbations where the fluctuations of both the bulk viscosity 
coefficient and the decay rate are invariant under infinitesimal diffcomorphisms. The gauge- 
invariant approach will be complemented with gauge-fixed descriptions that are particularly 
useful in connection with the intricacy of the numerical calculations. 

The present article is organized as follows. In Section 2 the basic framework of the 
analysis will be illustrated. Section 3 contains the gauge-invariant description of cosmo- 
logical perturbations in the presence both of a fluctuating decay rate and of a fluctuating 
bulk viscosity coefficient. In Section 4 the gauge-invariant evolution equations for the cur- 
vature perturbations will be introduced, while in Section 5 curvature perturbations will be 
discussed in the context of a gauge-fixed description. Section 6 contains various numerical 
examples clarifying the interplay between fluctuating bulk viscosity and fluctuating decay 
rate. Finally, Section 7 contains the summary of the main findings and the concluding 
discussions. For completeness, the analysis of the evolution of the vector modes has been 
included in an appendix. 
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2 Basic considerations 



Some physical quantities relevant to the problem at hand will now be introduced. As 
already pointed out, there are two independent physical effects that may compete in a 
mixture of relativistic fluids: the presence of a (fluctuating) decay rate and the presence of 
a (fluctuating) bulk viscosity coefficient. 

Consider a relativistic plasma formed by, at least, two species. The associated energy- 
momentum tensors can be written as 2 

^ = (Pa + Pa)«"Pa<r, 

rr = (Pb + PbK<-Pb«r- (2.i) 

If the fluids are decaying one into one another (for instance the a-fluid decays into the 
b-fluid), the covariant conservation equation only applies to the global relativistic plasma, 
while the energy-momentum tensors of the single species are not covariantly conserved and 
their specific form accounts for the transfer of energy between the a-fluid and the b-fluid: 

V M rr = -Tg l/a u a (p a + p a ), 

V.Tf = Tg va u a (p a + Pa ), (2.2) 

where the term T is the decay rate that can be both space- and time-dependent; in Eqs. 
(2.2) u a represents the (total) peculiar velocity field. Owing to the form of Eqs. (2.2), it 
is clear that the total energy-momentum tensor of the two fluids, i.e. T^ t = T£ v + T^ u is 
indeed covariantly conserved. 

Equations (2.2) can be easily generalized to the description of more complicated dynam- 
ical frameworks, where the relativistic mixture is characterized by more than two fluids. 
Consider the situation where the a-fluid decays as a — > b + c. Then, if a fraction / of the 
a-fluid decays into the b-fluid and a fraction (1 — /) into the c-fluid, Eqs. (2.2) can be 
generalized as 

v„rr = -rV Q u Q ( Pa + p a ), 
v„Tjf = /r<r*« a (p a + Pa ), 

v„ir = (i - f)r<r* u a ( Pa + Pa ), (2.3) 

and so on. 

2 The signature of the four-dimensional metric adopted here is mostly minus, i.e. (+,—,—,—). Latin 
subscripts are used to distinguish the pressures and energy densities of different fluids; tensor indices are 
instead denoted by lowercase Greek letters. 
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It is appropriate to recall here that the authors of Ref. [19] choose to parametrize the 
coupling between the fluids of the mixture through a more generic function. For instance, 
in the notations of [19], the source terms of Eqs. (2.2) will simply be Q v & and bound by 
the condition 3 Q v & + Q£ = 0. The latter condition follows from the covariant conservation 
of the total energy-momentum tensor of the mixture and it is also satisfied by the source 
terms in Eq. (2.2). The second condition satisfied by the source terms is that their spatial 
component has to be zero when evaluated on the Friedmann- Robertson- Walker background. 
In the present approach, the spatial part of the source terms given in Eq. (2.2) automatically 
vanish on the Friedmann- Robertson- Walker background since g^ u u^u u = 1 (implying Ho = 1 
and Ui = 0). 

The second relevant dissipative effect that can enter the energy-momentum tensor with- 
out spoiling the isotropy of the background geometry, as discussed in the previous section, 
is the one parametrized in terms of bulk viscosity. Following the notation of Eqs. (1.1) and 
(1.2), the effective energy-momentum tensors of the two fluids can be written as 

Hence, the Einstein equations can be written as 

K - = SttGT;, (2.5) 

where, 

= T£ V + T£ V . (2.6) 

According to Eqs. (1.1), (1.2) and (2.1), T£ v and will have an inviscid background 
value (characterized by and p a ,b) and a viscous correction parametrized, at the level of 
the background, by £ a and £ b . 

In the following the background geometry will be taken to be of Friedmann-Robertson- 
Walker (FRW) type. Moreover, it will be assumed that the background geometry is spatially 
flat so that the line element can be written as 

ds 2 = dt 2 - a 2 (t)dx 2 , (2.7) 



where t is the cosmic time coordinate. Therefore, Eq. (2.5) implies' 



i 



H 2 = —p, (2.8) 

H = -4ttG(p + V), (2.9) 
3 See, for instance, Eqs. (5.1), (5.2) and (5.3) of [19], page 28. 

4 We denote by an overdot a derivation with respect to the cosmic time-coordinate t. 
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where H = a/a is the Hubble rate; p = p a + pb and p = p a + Pb denote the total (inviscid) 
energy and pressure densities while V is the effective total pressure containing the viscous 
contribution, i.e. 

V = p - = (p a + p b ) - 3(£ a + ?„) H. (2.10) 
In the following sections, the total background viscosity coefficient will be denoted by £, i.e. 

e = ? a + e b + -, (2.H) 

where the ellipses stand for the possible contribution of other fluids. The definitions of the 
effective pressures for each fluid follows from Eq. (2.10) and it is 

V a = p a - 3? a H, V h = Ph - 3? b H. (2.12) 

Notice, furthermore, that £ (and, similarly £ a and £ b ) denotes the background values of 
the bulk viscosity coefficients. In different physical examples (see e.g. [4, 7] and also [11]) 
the bulk viscosity coefficient is proportional to a power of the energy density. If the a- and 
b-fluids would not be allowed to exchange energy and momentum, the mixture could be 
described by a single total fluid with effective pressure V . If the decay rate does not vanish, 
Eqs. (2.2) must be generalized so as to include both decay rate and viscous contribution: 

V^r = -Tg™u a (p a + p a ), (2.13) 
V M 7T = rg» a u a (p h + p h ). (2.14) 

In the case of a spatially flat FRW background, denoting with T the homogeneous compo- 
nent of the decay rate, Eqs. (2.13) and (2.14) reduce to 

Pa + 3fT(p a +Pa) + r(p a +p a )=0 (2.15) 

p h + 3H(p b + V h )-T(p a + Psi ) = 0. (2.16) 

Notice, once more, the difference between the effective pressures (i.e. V & and Vb defined in 
Eq. (2.12)) and the inviscid pressures (i.e. p a and pb). 

In the following sections we are going to study the evolution of the inhomogeneities 
in the class of background models introduced in the present section. Therefore, while the 
background quantities will be only function of the time coordinate, the fluctuations will be 
functions of both time and space. 

In the present paper the bulk viscous effects will be parametrized in the context of the 
Eckart approach [47] (which is the one also followed in [1]). It must be mentioned, that 
this approach is phenomenological in the sense that the bulk viscosity is not modeled on 
the basis of a suitable microscopic theory. For caveats concerning the Eckart approach see 
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[48] (see also [13] and references therein). This Eckart approach, however, fits with the 
phenomenological inclusion of a fluid decay rate that has been also considered recently for 
related applications to cosmological perturbation theory. Possible generalizations of the 
considerations reported in this paper to causal thermodynamics will not be examined here 
(see [40] for a covariant approach to this problem). 

3 Cosmological perturbations for imperfect fluids 

As is well known [29], density and metric perturbations may obey formally different evolu- 
tion equations, depending upon the specific gauge used to perform a given calculation. Once 
a specific gauge-invariant quantity is computed in a given gauge, it will be, by definition, 
the same in all the other coordinate systems. Not all coordinate systems are equivalent for 
specific applications, since one particular gauge may turn out to be more convenient for 
practical reasons. 

In the following the evolution equations of the fluctuations will be derived without 
choosing a specific gauge. This strategy has a twofold advantage. Firstly the full gauge- 
invariant evolution equations can be swiftly derived. Secondly, the most useful gauge choices 
can be selected by comparing the different forms of the equations in different coordinate 
systems. To begin with, the scalar modes of the Einstein equations will be first derived. We 
will the move to the evolution equations of the hydro dynamical quantities, i.e. the pressure 
and density fluctuations and the peculiar velocity of each fluid. While in the present and 
in the following sections the attention will be focused on the scalar modes of the geometry, 
some considerations on the evolution of the vector modes are collected in the appendix. 

3.1 Einstein system 

The scalar perturbation of Einstein equations can be formally written as 



where S s denotes the first-order scalar fluctuation of a given quantity and where Q£ denotes 
the Einstein tensor. Recalling Eq. (1.1) the scalar fluctuation of the effective energy- 
momentum tensor of the fluid sources then becomes, in general coordinates 




8nG5 s T> 



(3.1) 




= u v [{8p + 8p) - 5i V X 1I X - £ S(V x u x )] 
+ (5 Ufl u v + u» 5u»)[(p + p)- lV x u x ] , 



(3.2) 
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where the (total) velocity field satisfies g^u^Uy = 1 and where the bar means that the 
corresponding quantity has to be evaluated using the background geometry. To write Eq. 
(3.2) in explicit terms one has to make explicit the first-order fluctuation of the covariant 
derivative of the velocity field, i.e. 5(V\u x ). Recalling, in fact, that V\u x = d\u x + Tp a u a 
we have that 

5(V x u x ) = d x 5 u x + ST^vf* + T%Su a , (3.3) 

where Tp a and SF^ a are, respectively, the background and the perturbed Christoffel con- 
nections. 

Without performing any gauge choice, the scalar fluctuations of a conformally flat metric 
with line element 

ds 2 = g^dafdx" = a 2 (r) [dr 2 - dx 2 ] , (3.4) 

are parametrized by four independent functions so that the entries of the perturbed metric 
can be written as 

2 

5 s g 00 = 2a 2 4>, 5 s g 00 = ^ 



S s9ij = 2a 2 - didjE) , 5 s gV = - (# y - 8 l &> E) , 



5 s g 0i = -a 2 d % B, 5 s g 0i = -^B; (3.5) 



1 

note that Sij is nothing but the Kroneker's delta. 

In Eq. (3.5) the conformal time parametrization of the background geometry has been 
assumed and the evolution equations of the background can be written as 5 

n 2 = ^fa 2 P , (3.6) 

7i 2 -H' = 4irGa 2 (p + V), (3.7) 

p' a + 3H(p a + V a ) + aT{ Pi ,+p a ) =0, (3.8) 

p' h + 3H( Ph + Vb) ~ aT(p a + p a ) = 0, (3.9) 

where H = a' /a. In Eqs. (3.7)-(3.9) the definition of the total and partial shifted pressures 
follows from Eqs. (2.10) and (2.12). Since g^u^Uy = 1, it also follows, from Eqs. (3.4) and 
(3.5) that 

Ho = a, 5u° = —<j)/a. (3.10) 

Equations (3.3) and (3.5) lead to the following explicit expressions of the fluctuations 
of the effective energy-momentum tensor 

5%° = Sp, (3.11) 



'We denote by a prime a derivation with respect to the conformal time coordinate r. 
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67q = [ P + P -£h)v\ (3.12) 



a 



5T-> = -8{{5p-?^5t-^[0-Z(il)' + H<j>)+V 2 E'\y (3.13) 

where the following notation for the peculiar velocity field and for its divergence has been 
chosen: 

uoSu* = v\ dy = 9. (3.14) 

Since no specific gauge foxing has been invoked in Eqs. (3.11)-(3.13), Eq. (3.1) allows the 
determination of the Hamiltonian constraint 

VV - HV 2 {B - E') - 3H{tp' + H<f>) = 4irGa 2 5p, (3.15) 

and the momentum constraint 

V 2 (V>' + H<f>) + (H 2 - H')(V 2 B + 0) = 0. (3.16) 

In the perturbed geometry (3.5), the (i = j) and (i / j) components of Eq. (3.1) 
become, respectively 

i>" + (TC 2 + 2H')<f> + n{4>' + + ^V 2 [(<P -il>) + (B- E')' + 2H{B - E')] 

= 4irGa 2 ^5p - 3^5£ - ^[9 - 3(Y>' + H</>) + V 2 E']j, (3.17) 

-^di&^E' - B)' + 2H(E' -B) + (i;- 0)] = 0. (3.18) 

In Eq. (3.18) the possible contribution of the anisotropic stress has been neglected. 

For infinitesimal coordinate transformations, the fluctuations of both the metric and 
the effective energy-momentum tensor change. Hence, it is always convenient to define 
appropriate gauge-invariant quantities. A possible set of gauge-invariant quantities is then 

<S> = ( f )+ ( y B-E'y + H(B-E') (3.19) 

* = -0 - H(B - E') (3.20) 

for the metric and 

5 Pg = 5p + p'{B - E') (3.21) 

5 Pg = 5p + p'(B -E') (3.22) 

= 9 + V 2 E' (3.23) 

for the (total) fluid sources. In Eqs. (3.21) and (3.22) the subscript "g" means that the 
corresponding quantity is invariant under infinitesimal coordinate transformations, i.e., for 
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short, gauge-invariant. This notation will also be employed, when needed, in the following 
sections. Equations (3.19) and (3.20) are the gauge-invariant fluctuations of the metric 
sometimes called Bardeen potentials [29]. From Eqs. (3.19) and (3.20), it is clear that in 
the gauge E = and B = <!> and ^ coincide with the longitudinal fluctuations of the 
metric, i.e. (ft and ift [30]. 

Since the bulk viscosity coefficient is general time dependent, its gauge-invariant fluctu- 
ation is defined as 

Z = 5Z+t(B-E'). (3.24) 

Following the usual procedure, the various fluctuations both of the geometry and of the 
sources, i.e. ift, (ft, bp , 5£ and so on, can be expressed in terms of their gauge-invariant 
counterparts, i.e. Sp g , E and so on. These expressions can then be inserted into Eqs. 

(3. 15), (3. 16) and (3. 17), (3. 18). The resulting set of equations will now contain only gauge- 
invariant quantities. The full system will then be formed by the gauge-invariant counterpart 
of Eqs. (3.15) and (3.16): 

V 2 f - 3W(*' + W$) = 4irGa 2 5 Pg , (3.25) 
V 2 (*' + H<S>) + (H 2 -H')e = 0, (3.26) 

and by the gauge-invariant counterpart of Eq. (3.17): 

+ + 2*') + (7i 2 + 2W')$ = 47rGa 2 |<5pg - ™E - J[0 - 3(*' + (3.27) 
Equation (3.18) implies that <3? = 'I' in the absence of anisotropic stress. 

3.2 Interacting viscous mixtures 

The perturbed form of the covariant conservation equations can be written as 

<5(V^) = - [ST g ua u a (p a + p a ) + r 5g va u a (p a + p & ) 

+ Tg ua Su a (p a + P,) + rg ua u a (5 Pll + 5p a )}, (3.28) 
«5(V M T ( ^) = [5Tr a u a (v ! , + p a )+T5g va u a {p & + p & ) 

+ rr Q <WPa + Pa)+rr Q M^a + ^a)}. (3-29) 

In [19] the authors choose to parametrize the rate of momentum transfer (in the pertur- 
bation equations) through a more generic functions whose properties follow, in the present 
approach, from the form of the coupling already discussed after Eqs. (2.2) (see, for instance, 
Eq. (5.18) of [19]). 
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From Eq. (3.28), without fixing the gauge, we obtain the following general expression 
from the (0) component 

5p a + m (5 Pa + 5 P& ) - 9^ 5£a + (p a + p a ) 9 a - 3 (p a + P a ) $ 

-3 ^ W[0 a - 3(V' + W$ + V 2 ^'] + ( Pa + P a ) V 2 £' 

a 

= -a[(5r ( Pa + Pa ) + r( Pa + p a ) + r (5 Pa + <5pa)] , (3.30) 

while from Eq. (3.29), still for the (0) component, the result will be 

H 2 

5p' h + m {5 Ph + <5p b ) - 9— 5C b + (p b + P b ) h - 3 (p b + P b ) ^' 
-3 ^ H [9 h - 3 (V' + H<j>) + V 2 E'} + (p b + P b ) V 2 £' 

= a[5T (p a + Pa) + 4> r (p a + Pa) + T (5p a + <5p a )] . (3.31) 

To derive Eqs. (3.30) and (3.31) from Eqs. (3.28) and (3.29), Eqs. (3.8) and (3.9) have 
been repeatedly used. 

The (i) component of Eqs. (3.28) and (3.29) leads, respectively, to the evolution of the 
peculiar velocity of the a-fluid: 



+ 



m + 



p', + K 



Pa+Pa 



+ V 2 (B' + <t>) + 



V 

n+ a 



Pa + Pa 



o „ V 2 (5p a 
Pa + Pa 



3HV 2 5^ £ a [V 2 ^ a - 3V 2 (^' + H<f>) + V 4 £'] = _ af P, + Pa 



a (p a + Pa) a(Pa + Pa) 

and to the evolution of the peculiar velocity of the b-fluid: 



+ 



m + 



Pb + Pb 



e h + v 2 {B' + 0) + 



n+ b 



Pb + Pb 



Pa + Pa ' 



„2 „ V 2 5 Ph 
X7 2 B + 



(3.32) 



Pb + Pb 



3WV 2 <5£b £ b [V 2 ^ b - 3V 2 (V + H$) + V 4 £'] _ ar p a +Pa g (3 33) 



a (Pb + Pb) 



a(Pb + V h ) 



Pb + V h 



Equations (3.30), (3.31) and (3.32), (3.33) can easily be written only in terms of quantities 
that are explicitly gauge-invariant. To this purpose it is useful to point out that the gauge- 
invariant fluctuations related to the bulk viscosity coefficients and to the decay rate can be 
written as 



S (a,b) = <5£(a,b) + C(a,b)(- B ~ E ')i 

5T g = 5T + T'(B - E'), 



(3.34) 
(3.35) 



where the subscripts in Eq. (3.34) imply that the gauge-invariant fluctuation has the same 
form for both the a- and the b-fluid. 
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Using Eqs. (3.34) and (3.35), Eqs. (3.30), (3.31) and (3.32), (3.33) can be written in 
fully gauge-invariant terms. More specifically, Eqs. (3.30) and (3.31) lead, respectively, to 

5p' ga + 3H(5 Pga + 5 Pga ) - 9—E a + (p a + P a )0 a - 3(p a + P a )*' 

-^ a [e a -3(^ + = 



and to 



-a[ST g (p a + p a ) + r((5p ga + 5p ga ) + T(p a + p a )$] 



<5pl b + 3W(<5p gb + 5 Pgh ) - 9—E h + ( Ph + P b )0 b - 3(p b + V h )V 



(3.36) 



— W£ b [e b -3(tf' + W$)] = 
a 

a[<5r g (p a + p a ) + r((5pg a + <5p ga ) + r(p a + p a )$], 

while Eqs. (3.32) and (3.33) become 

V 2 <5p g 



(3.37) 



©a + 



+V 

p a + V a 



e a + 



a(p a + V a 



^ _e a [v 2 e a -3v 2 (M/' + Hci,)] 

" v ^ a 



ar Pa + Pa @ 
Pa + Pa 



e' b + 

-3- 



m + 



Pb + V h 



-V 2 H b 



a(p a + Pa) 

Pb + T> h 

| b [V 2 B b - 3V 2 (^ + Tin _ ^ Pa + Pa 



ar- 



e. 



(3.38) 



(3.39) 



a(Pb + Vb) a(Pa + Vb) " Pb + V h 

In Eqs. (3.36), (3.37) and (3.38), (3.39), 5p ag and <5p bg denote the gauge-invariant energy 
density fluctuations of the two species. A similar notation has been employed for the fluctu- 
ations in the pressure densities (i.e. 5p ag and <5p bg ) and for the gauge-invariant fluctuations 
of the velocity field (i.e. a and 0b)- A useful remark is that Eqs. (3.38) and (3.39) imply 
a transfer of momentum between the a-fluid and the b-fluid. 

The total gauge-invariant fluctuations are easily related to the gauge-invariant fluctua- 
tions of the single species, i.e. 



5p g = 5p ga + 5p gh , 
5 Pg = 5p ga + 5p gh , 

(p + p)0 = (p a + p a )0 a + (p b + Pb)©b, 

ee = e a ©a + ? b e b , 

E = E a + E h . 



(3.40) 
(3.41) 
(3.42) 
(3.43) 
(3.44) 



From Eqs. (3.42) and (3.43), it follows that the total velocity field can also be written as 

(p + V)@ = (p a + P a )0 a + (p b + P b )e b , (3.45) 
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where, we recall, V = V & + V\,. 

It is clear that summing up Eq. (3.36) and Eq. (3.37), the temporal component of the 
covariant conservation equation for the total fluid is recovered, i.e. 

5p' g + 3H(5 Pg + 5 Pg ) - 9^ ~ + (p a + P a )G 

-3(p + P)*' - -W?[6 - 3(V + = 0. (3.46) 

Similarly, summing up Eq. (3.38) and Eq. (3.39), the spatial component of the covariant 
conservation equation for the total fluid becomes 



e' + 



m + 



(p' + V) 



p + v 



e + — £| + v 2 $ 



a(p + P) p + V v ; 

As a final comment, it could be noted that while, in the case of the background, the effect 
of the bulk viscosity coefficient is to renormalize the effective pressure densities, at the 
level of the fluctuations this simplification is no longer possible as is evident, for instance, 
from Eq. (3.27). In a complementary perspective, the presence of spatial fluctuations in 
the bulk viscosity coefficient induces some sort of intrinsic non-adiabatic pressure density 
variation. Each of the fluids of the mixture will have an inviscid contribution and a viscous 
contribution. In the limits £ ab —> and H aj b — > 0, the fluid becomes inviscid both at the 
level of the background and at that of the fluctuations. The inviscid contribution can be 
characterized, in the simplest case, by two barotropic indices and by the two sound speeds, 
i.e. 

w a , h - - — , c sa - — — , c sb - — — , (3.48) 

Pa,b <jPga OPgh 

where, in the first equation, the subscript means that the relation holds independently for 
each fluid. We are now in a condition to derive the evolution equations of the curvature 
and entropy perturbations in the case when both the decay rate and the bulk viscosity 
coefficient experience some spatial variation. 



4 Curvature and entropy perturbations 

The curvature fluctuations usually discussed in the theory of cosmological perturbations 
[29, 49, 50, 51] are indeed perturbations of the spatial curvature on comoving orthogo- 
nal hypersurfaces (see, for instance, [52]). Gauge- invariant curvature fluctuations can be 
defined, in terms of the Bardeen potentials, as 

= -(*-W g ), (4.1) 



n = - 



TL 



n 2 -n 



7 (w$ + ^') 
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where, in the second equality, we used the notation 6 6 = V 2 U g . Note that the second 
equality follows by using the momentum constraint. Suppose, indeed, to be in the comoving 
orthogonal gauge where, by definition (vc + -Be) = 0. In this gauge, the fluctuations of the 
spatial curvature are related to ipe- Suppose now to shift by an infinitesimal quantity eo 
the conformal time coordinate tq, i-e. 

re -»■ t = r c + e . (4.2) 

Under the coordinate shift (4.2), both V'c an d (i>c + ^c) transform as 

(v c + B c ) ->■ v + B = v c + B c + eo, (4.3) 
ipe i> = i>G + We . (4.4) 

Consequently, the expression of tpc in the new coordinate system defined by Eq. (4.2) is 
simply 

if> c = t p-n(v + B) = (y-HV g ). (4.5) 

In the first equality of Eq. (4.5) the value of eo is fixed by Eq. (4.3). The second equality 
of Eq. (4.5) follows, instead, from the definition of gauge-invariant fluctuations of Eqs. 
(3.20) and (3.23). Clearly, up to a sign (this is a matter of conventions) Eq. (4.5) is exactly 
Eq. (4.1). While in the comoving orthogonal gauge, 1Z is related to fluctuations of the 
spatial curvature, in a different coordinate system 1Z will have the same numerical value 
(by gauge-invariance) but will not necessarily be related to curvature fluctuations. 

There is also another physical quantity that plays an important role in our considera- 
tions, namely the curvature fluctuation on uniform density hypersurfaces (see, for instance, 
[52]). The fluctuations of the (total) curvature on uniform density hypersurfaces have a 
gauge-invariant interpretation and they can be defined as 



C = -(* + H-^| (4.6) 

As in the case discussed in Eq. (4.1), ( coincides with a perturbation of the spatial curvature 
in the uniform density gauge (see, for instance, [52]). In a diffferent coordinate system £ will 
keep the same numerical value, but its physical interpretation will be different. Suppose 
to be in the coordinate system defined by the condition 5/?d- This is the so-called uniform 
density gauge. Consider then curvature fluctuations on the uniform density hypersurfaces 
V>d- Under an infinitesimal shift of the time coordinate, ips an d 5pn change as 

V'D > ip = ipD + He , (4.7) 
5pB -> 5p = 5pr, - p'e . (4.8) 



6 In the following paragraphs of the present section the notation 8 — V 2 u and = V 2 V g will be employed; 
these definitions follow naturally from the conventions of Eqs. (3.14) and (3.23). 
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Equations (4.7) and (4.8) imply that, in the new coordinate system, 

V, D = V + W^ = * + W^, (4.9) 
P P 

where the first equality follows from Eq. (4.8) and the second from Eqs. (3.20) and (3.21). 
Again, up to a sign, Eq. (4.9) coincides with Eq. (4.6). This means that Q parametrizes 
fluctuations in the spatial curvature computed on uniform density hypersurfaces. In a 
different gauge the numerical value of Q will remain the same but its interpretation may be 
completely different. 

By taking the difference of Eqs. (4.1) and (4.6) and by using the Hamiltonian constraint 
in its gauge-invariant form, i.e. Eq. (3.25), the following relation between TZ and £ can be 
obtained: 

n '<- i2S( P + vy (410) 

where, as usual, V is the total shifted pressure of the mixture. Equation (4.10) implies that, 
up to Laplacians of ^, TZ ~ Q. Equation (4.10) generalizes the relation between TZ and ( 
to the case when the bulk viscosity coefficient is non- vanishing. In fact, in the denominator 
at the right-hand side of Eq. (4.10) there appears V, i.e. the sum of inviscid and viscous 
contributions to the background pressure. 

As elaborated, Eq. (4.6) describes curvature fluctuations on uniform density hyper- 
surfaces. In the same spirit, it is also possible to define the curvature fluctuations on the 
hypersurfaces where the energy densities of either the a-fluid or the b-fluid are uniform, 
namely: 

Ca = -(*+W^), (4.11) 

Ch = -U + H *£pL\ (4.12) 
V Pb / 

where, according to Eqs. (3.36) and (3.37), 5p ga and Sp g b are the gauge-invariant density 
fluctuations of each fluid of the mixture. In equivalent terms, Eqs. (4.11) and (4.12) can 
also be interpreted as related to the density fluctuations on the hypersurfaces where the 
curvature is homogeneous (this point of view will be particularly useful for the considerations 
of the following section) . 

From the gauge-invariant definitions given in Eqs. (4.6) and (4.11), (4.12) it is easy to 
show that / f 

C=^Ca + ^Cb- (4.13) 

While the weighted sum of £ a and Cb is related to (, the difference between them is related to 
the entropy fluctuations. More specifically, the relative fluctuations in the specific entropy 
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? can be written as 

5 = ^ = -3(Ca-Cb). (4.14) 

In the limit of vanishing decay rate and vanishing bulk viscous stresses (for two perfect, 
inviscid, non- interacting fluids), Eq. (4.14) becomes 

s = orky " (i+k) • (415) 

where S ga = 5p ga /p a and 5 g b = Sp g b/Pb are the gauge-invariant density contrasts; w a and 
Wb are the barotropic indices for the two fluids of the mixture. Equation (4.15) coincides 
with the definition of Ref. [19] and it applies, for instance, in the discussion of CDM- 
radiation isocurvature mode in the absence of decay rates. In more general situations, for 
instance when decay rates are included, Eq. (4.14) does not reduce to Eq. (4.15) [27]. 
Recalling Eqs. (4.11) and (4.12), Eq. (4.14) can also be written as 

s = m (^_Spf>\ ( 4.16) 

Let us now define the quantities relevant to the evolution of curvature perturbations. 
The total pressure density appearing in Eq. (3.27) can always be split as 



where the two subscripts imply that the two relative variations at the right-hand side 
should be taken, respectively, at constant entropy and energy densities. Thus, to perform 
the variation at constant (total) energy density means that 5p ga = —Sp g b- Similarly, to 
perform the variation at constant q means that 5$ = 0, i.e. from Eq. (4.16) 

<fy>ga _ 5p g b 



Pi Pb 



(4.18) 



Notice also that, for practical purposes, it is often useful to denote the whole second term 
at the right-hand side of Eq. (4.17) as 5p n ad- 

In the case of two fluids the total speed of sound and the non-adiabatic pressure density 
variation are: 

n 2_(S_Pg\ _ CsaPa + jjX (a i q\ 

s ~Up g )r p' a +p' h ' { } 

* nad " {^Jf" ~ ~ H(pi + P > h ) (Ca " Cb) ' (420) 

where the speeds of sound in the two fluids of the mixture have been explicitly introduced. 
Recalling the connection between ( and the weighted sum of Ca and Cb> i-e- Eq. (4.13), it 
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is also possible to write 5p na d in an alternative useful form: 

froad= (C ' b ~^ a)p/a (Ca-C), (4-21) 

where, according to Eq. (3.48) the speeds of sound refer to the inviscid contribution to the 
total energy-momentum tensor. Thus c 2 b = w b and c 2 a = w a . 

Taking now the difference between Eqs. (3.27) and (3.25), and recalling Eqs. (4.1) and 
(4.17), the evolution of 1Z can be obtained easily, and it turns out to be 

_,, 3Tt t n Ti. r 37i 2 „ 

a(p + V) p + V a(p + V) 

+^^—0 ^ V 2 *. (4.22) 

We can also write the above expression directly in terms of £, i.e. as the density contrast 
on uniform density hyper surf aces. Recalling Eq. (4.10) we can also write 

n> = t> ^ n(^ + i) v2 „ 

^ 12irGa 2 (p + V) 12vrGa 2 (p + V) 

n ^ v 2 ^ + ^ -m v2 ^- ( 4 - 23 ) 



4TrGa 2 (p + V) 2 a a(p + V) 

by simply taking the conformal-time derivatives of the left- and of the right-hand sides. 
Inserting Eq. (4.23) into Eq. (4.22) and using the gauge-invariant form of the momentum 
constraint given in Eq. (3.26), we simply obtain 

37Y — / > Ti, . 3Ti , n ,_ — / T . 

a(/9 + 7 3 ) p + a(p + V) 



+G 



£W 1 



a(p + 7>) 3 



-V 2 ($ - *) - — -i -V 2 ^. (4.24) 



l2^Ga 2 {p + V) v 7 a(p + P) 



Concerning Eq. (4.24), a few comments are in order: 

• the first three terms on the right-hand side of Eq. (4.24) dominate over length scales 
larger than the Hubble radius; the remaining terms in the second line of Eq. (4.24) 
are sub leading in this limit; 

• when £ — > and S — > (i.e. in the case of vanishing homogeneous and inhomogeneous 
viscosity), Eq. (4.24) reduces to the usual form of the evolution equation for £ (note, 
in fact, that in this limit V — > p); 

• since <5p na d contains both £ and £ a (being proportional to the fluctuations of the 
entropy density), the evolution of £ is completely specified only when coupled to the 
evolution of £ a . 
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Therefore, to have a self-consistent system of differential equations for length scales 
larger than the Hubble radius, it is mandatory to derive an evolution equation for £ a . This 
can be deduced by inserting Sp ga (obtained from Eq. (4.11)) into Eq. (3.36). The result of 
this algebraic operation can then be written as 



(Ca + *)' + (Ca + *) 

+ , ^a , 

a Pa Pa 



[-( 









+ 3H(c sa + l) + ar(c sa + l) 



a Pi 



= a(p a + p a )r— 

Pa, 



H'-H 2 

h 2 



K 



H 

_ „ (Pa + Pa) ^ 
1 

Pa 



H 2 - Ti! 

n 



where 



(4.25) 
(4.26) 



is the fractional fluctuation of the decay rate. 

Equation (4.25) can be simplified by using the evolution equations of the background 
and, in particular, the background conservation equation for the a-fluid, i.e. Eq. (3.8). 
This observation allows us to eliminate the various terms proportional to \I> / and to \£. The 
second achievable simplification is to trade 1Z for (. This operation, according to Eq. (4.10), 
will bring some extra terms proportional to V 2 ^. In the case where the inviscid part of the 
a-fluid has a barotropic index w a , the speed of sound is given by c 2 a = w a . In this case the 
quantity 

H(l + w a )p a 



.9a 



Pa 



(4.27) 



can be defined. 

Taking into account all these simplifications, Eq. (4.25) can be written as 



Ca + 



9H 2 £ a 



5a 



m 2 - 

OPa 



L a p a 

e a + 3 



OPa 



5a 

n 2 -n' 

H 



+ 3 Z% V 2 * 

apk 
H'-H 2 



H l 



(4.28) 



In Eq. (4.28), several terms are subleading for length scales larger than the Hubble radius. 
Thus, we can rewrite Eqs. (4.24) and (4.28) in the cosmic time coordinate, keeping only 
the dominant terms for length scales larger than the Hubble radius. The result is given by 

h , m 



(P + 'P) p + v 



(p + r) 



(HE + £*) 



(4.29) 
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and by 



Ca + 



g^2^a _ 5a 
Pa 9a 



(4.30) 



Equations (4.29) and (4.30) can be used, for instance, to study numerically the evolution 
of fluctuations whose wavelength is larger than the Hubble radius, with two caveats. The 
first caveat is that g a has to be always non-singular. This may not always be the case when 
a viscous contribution is present. In this case, instead of studying the evolution of Ca it 
will be more practical to deal directly with <5 ga , i.e. the gauge-invariant density contrast. 
The second remark is that, apparently, Eqs. (4.29) and (4.30) require the independent 
determination of The evolution equation for \I> can be directly obtained in terms of ( 
from Eqs. (4.1) and (4.10). As will be shown in the next section, the coupled evolution of Ca 
and ( can also be easily integrated directly in a specific gauge, i.e. the so-called off-diagonal 
gauge [53] (see also [54] and [52]). 

It is appropriate to remark here that in the case of homogeneous decay rate (i.e. 5r g = 0) 
and in the absence of homogeneous and inhomogeneous viscosity (i.e. £ a = £ b = and 
H a = Sb = 0), Eqs. (4.29) and (4.30) become simply, in units 8ttG = 1, 

C= : MOCa), ( 4 - 31 ) 

P + P 

Ca = r 5a -^(Ca-C)- (4-32) 

Equation (4.31) can be deduced from Eq. (4.29) using Eq. (4.20) in the specific case 
Cg a = w a and c^ h = w^. To get to Eq. (4.31) it is also useful to recall that £b can be 
always eliminated through Eq. (4.13). Note that in Eq. (4.31), at the right-hand side, the 
denominator reads (p + p) and not (p + V) since V — > p in the limit £ = 0. Equation (4.32) 
can be obtained from Eq. (4.30) by making explicit the expressions depending on g a and 
by using Eq. (2.15) always in the case a = w a . 



5 Off-diagonal gauge description 

In the off-diagonal gauge, the diagonal entries of the perturbed metric are set to zero, i.e. 
5 s gij = in Eq. (3.5). More precisely, in the off-diagonal gauge, 

V'od = 0, E od = 0. (5.1) 
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From a generic perturbed line element where all the entries of the perturbed metric do 
not vanish, the off-diagonal coordinate system can be reached by performing a coordinate 
transformation 



T -> f = r + e , 

Xi — »■ Xj = Xj + <9«A. (5.2) 
The shift (5.2) induces a transformation in the S s gij components of the perturbed metric: 

fp -»■ ^od = ip + 1-Ce , 

E^E od = E- A. (5.3) 

The gauge parameters are completely fixed to A = i£ and eo = —tp/H since, as introduced 
above, V'od = -^od = 0. In the off-diagonal gauge the curvature is then uniform, i.e. ip Q d = 0. 
This is the reason why the off-diagonal gauge is also correctly named uniform-curvature 
gauge. 

In the off-diagonal coordinate system, the gauge-invariant variables Ca and Cb defined 
in the previous section are simply proportional to the relative density contrasts of the two 
fluids, i.e. 

Pa Pa 

C b = -H 5 ^ = ~H^5 h , (5.4) 
Pb Pb 

where <5 a and 5b are the density contrasts computed in the off-diagonal gauge. 

The evolution equations in the off-diagonal gauge can be obtained directly from the Eqs. 
(3. 15), (3.16) and (3. 17)-, (3.18) by setting E and ip to 0. With this strategy in mind the 
Hamiltonian and momentum constraints read, respectively 7 , 

-HV 2 B - m 2 4> = 4irGa 2 5p, (5.5) 
HV 2 cf> + (H 2 - H')V 2 B = -4irGa 2 (p + p)9, (5.6) 

while the (i = j) and (i / j) components of the perturbed equations of motion are, respec- 
tively, 

(H 2 + 2H')4> + H<f>' = 4TrGa 2 5p, (5.7) 
B' + 2HB + <f> = 0. (5.8) 



7 In the following, without ambiguity, <j>, B, 8 and all the other perturbation variables are meant to be 
evaluated in the off-diagonal gauge. 
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The covariant conservation of the perturbed energy-momentum tensor implies, instead, 
<5p' a + m{5p a + 5p a ) + (p a + Pa )9 a = -o(p a + Pa )T (s r + + S P- + Sp A , (5.9) 

V Pa + Pa, / 

5p' h + 3H(5 Ph + tfpb) + (p b + Ph )9 h = a(p a + p a )r(V + </> + ^ a j^M . (5.10) 

V Pa+Pa / 

In the off-diagonal gauge the explicit form of the evolution equations for £ and £ a (ampu- 
tated to the leading terms that are relevant for length scales larger than the Hubble radius) 
are: 

C = (*£a + ^b) + (fa + fb)C] + i(»b " ^a)Pa(Ca " C), (5-11) 

where we traded (p + V) for —2H, as follows from Eq. (2.9) (when units 8irG = 1 are used 
8 ). The evolution equation for ( a is then 

C a + (V^ - is) c a + — {HHu - h~u) = -r 5a (sr + ^c) ■ (5.12) 

V Pa 5a/ Pa V H z J 

To derive Eqs. (5.11) and (5.12) directly in the off-diagonal gauge, it is useful to notice 
that, in this coordinate system, the relation of £ to is particularly simple and it is given 
by 

T-/ 2 H 2 

< = -W=k* s T*- (5 - 13) 

As previously discussed, it is useful for some applications to study curvature fluctuations 
directly in the £, 5 a parametrization where the evolution equations read: 

C = ~^[H(SU + <%) + (!a + fb)C] - ^K" w a )(H Pa 5 a + p a C), (5.14) 
In In 

k + (?A " <5£a + la^c) = "^(1 + W a ) (V + ■ (5-15) 

Notice that Eqs. (5.11) and (5.12) can be swiftly deduced from Eqs. (4.29) and (4.30), 
by expressing the gauge-invariant fluctuations directly in the off-diagonal gauge. The only 
non-trivial expressions are 

H~+ty = n5£, (5.16) 



and 



HE a + = (5.17) 



The meaning of Eqs. (5.16) and (5.17) is that the two (gauge-invariant) quantities at the 
left-hand sides translate, in the off-diagonal gauge, simply into the spatial variations of the 
bulk viscosity coefficient. The right-hand sides of these equations are deduced from the 



In the following, without loss of generality, we will set 8nG = 1. 



22 



left-hand sides by simply using Eqs. (3.20)-(3.24) (or Eqs. (3.20)-(3.34)) and by setting, 
in the obtained expression, E = and tp = 0. 

For some applications, it is interesting to consider situations where H — > 0, at least for 
an interval of cosmic time. This may happen if the background solutions have, in some limit, 
a quasi-de Sitter evolution where H^ 1 may become very large. In this case the evolution 
equations in terms of ( is partially invalidated. The evolution equations in terms of (f> are, 
instead, fully valid and well defined in all their limits. 

From Eqs. (5.7), (5.8) and (5.9), neglecting the spatial gradients and using the cosmic 
time coordinate the relevant equations become 

4> + [zH + 2|^ = ^{5p + 3H[g a + Z h )cf> - (<5£ a + «$&)]}, 

9H 2 - - - 

^ + ~ + L<P] = -I> a + l)[Sr + 0], 

Pa 

p + H(3 + (f> = 0, (5.18) 

where aB = (5. 



6 Physical applications 

6.1 The case of vanishing background viscosity 

Two complementary situations may arise. In the first case the homogeneous component 
of the bulk viscosity coefficient vanishes (i.e. £ a = £ b = 0) but the inhomogeneous part 
does not vanish (i.e. either 5^ / or 5^ / 0). In this case it is plausible, from either 
Eqs. (4.29), (4.30) or (5.16), (5.17), that the evolution of curvature perturbations will be 
supplemented by a source term proportional to the spatial variations of the bulk viscosity 
coefficient. 

This situation is, to some extent, a bit arbitrary, but it closely analogous to the situation 
encountered in the case of a decay rate with spatial fluctuations, i.e. the toy model when 
the decay rate is allowed to be inhomogeneous but without a specific microscopic model 
supporting this interpretation [41]. In spite of this important caveat let us investigate first 
this (rather extreme) situation. 

Suppose, for simplicity, that £ a = £ b = but 8^ / 0. In this case, it can be shown 
that, in spite of the value of 5r, C increases linearly in time. In Fig. 1 the integration of 
Eqs. (5.11) and (5.12) is illustrated. In the plot at the left-hand side, the evolution of the 
energy densities of the background is reported for the case where the decay products are 
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Figure 1: The evolution of the background (left-hand-side plot) and of the fluctuations 
(right-hand-side plot) are illustrated in various cases of physical interest. The parameters 
are fixed in so that T/Hq = 10~ 3 , 5r = 1CT 7 , c5£b = 1CT 10 , &; a = 0. In this figure (and in 
the following) Hq denotes the value of the Hubble parameter at the initial integration time. 



given by a radiation fluid (i.e. = 1/3) and the a-fluid is given by pressureless matter 
(w a = 0). Similar plots can be obtained when w a ^ (for instance w a = 1 or w a = 1/10). 
In the present context it will often be interesting to fix to 1/3 and to leave w a free to 
vary form to 1. As is clear from Fig. 1, for t > T , the energy density of the a-fluid is 
exponentially damped thanks to the finite value of the decay. For curvature scales much 
smaller than the decay rate the radiation background dominates. In Fig. 1 (plot at the 
right-hand side) the evolution of £ is compared with the analytical approximation valid for 
large times. 

To get an analytical approximation for the evolution of £ let us start from Eqs. (5.11) 
and (5.12). We are interested in the case where £ a = £ b = 0. Moreover, as assumed above, 
we may also posit, for simplicity, w a = (which is the case numerically treated in Fig. 1). 
Under these assumptions, Eqs. (5.11) and (5.12) become, respectively 

C = ^Pa(Ca - C) " ^6t, (6.1) 

Ca = r-^ 5 a(Ca - C) - r 9sL 5 r . (6.2) 

To proceed further, notice that, from Eq. (4.27) we have, in the case w a = 

Hp a H 

5a = "77 = 3H + T (6 - 3) 
where the second equality follows from the first by explicit use of the (background) continuity 
equation for the a-fluid in the case w a = 0: 

p a + (3H + T)p a = 0. (6.4) 
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Since T is constant, from Eq. (6.3) it is clear that, for large times, 

5a ^ H/T. (6.5) 

This asymptotic expression compares very well with the numerical result. In Fig. 1 (plot 
at the left), the first thick line from the top corresponds to the numerical evolution of logg a 
while the (first) thin line (from the top) is the result obtained by approximating g a ~ H/T. 

Notice, furthermore, that for sufficiently large times (i.e. t > T 1 ) the evolution of H is 
simply given by 

H ~ — — -. (6.6) 

3(w h + l)t K ' 

This is because the contribution of the a-fluid to the Hubble rate vanishes exponentially for 

t > T 1 . In fact, from Eq. (6.4) direct integration leads to 

Pa(t)=Pa(to)(^)V^-*o). (6.7) 

Equations (6.4)-(6.5) and (6.6)-(6.7) allow to estimate analytically the time-dependent 
coefficients appearing in eqs. (6.1) and (6.2). In particular, by looking at Eq. (6.1), the 
coefficient of the first term at the right hand side is exponentially suppressed, for t > T 1 
since it contains p & (see Eq. (6.7)). On the contrary, the coefficient of the second term at 
the right hand side of Eq. (6.1) is finite and it is given by (see Eq. (6.6)) 



3H 



2 



2 H w h + l 

Since <S£ is constant in time, direct integration of Eq. (6.1) leads to 



(6i 



Wh + 1 \h 

where t±>T 1 . By keeping the dominant terms in the limit t 3> T Eq. (6.2) allows also 
to determine the evolution of Ca that turns out to be of the same order of (. Equation (6.9) 
is illustrated in the right-hand plot of Fig. 1 and it compares very well with the numerical 
solution. 



6.2 The case of non vanishing background viscosity 



The following example we ought to discuss is the one where £ a = e^fpl and = 0. In this 

case, 



5£a 



fya 
: 2^Pa 



L = e 



Pa 
2^ 



(6.10) 



25 



Moreover, in the off-diagonal gauge, the following chain of equality holds 

s 2 2H^/pT v ; 

It will be relevant, to physical applications, to analyse the case when the a-fluid decays 
into radiation (i.e. = 1/3) with different values of the homogeneous decay rate and 
different initial values of w a . When the a-fluid decays, for r ~ H, energy is transferred to 
the b-fluid. Different values of e affect the overall normalization of the energy density of the 
a-fluid. However, since the a-fluid decays exponentially and the b-fluid decays as a~ 3 ( Wb+1 ^ 
(where a is the scale factor) the a-fluid will quickly be subdominant. Different values of 
the initial decay rate (i.e. T/Hq < 1) may shift the onset of the b-dominated phase. It is 
relevant to stress, at this point, that initial conditions for the evolution of the background 
are set in such a way that pb(^o) = (where to is the initial integration time). Initially, 
around to, the a-fluid dominates and it can easily be shown, recalling that £ a = e^fp~i that 
the approximate solution for the evolution of the background, in this regime, is given by 

H ®~w i? o 7Tt7 ' Pa (t)~3H 2 (t), t <t<T-\ (6.12) 

[3(u; a + 1) - 3V3eJt 

In fact, for t < T 1 the a-fluid dominates the energy-momentum tensor, therefore from Eqs. 
(2.8)-(2.9) and (2.15)-(2.16) (written in units 8irG = 1 that are the units employed in our 
numerical codes) we can obtain the following equation 

H = -^[{w a + l)-V3e]H 2 (6.13) 

whose integral with respect to cosmic time is exactly Eq. (6.12). For to<t < T 1 , the scale 

factor a(t) expands (i.e. d > 0) in a decelerated way (i.e. a < 0) and H never goes to zero 

iff e < (w a + l)/\/3. Notice that, in the limit e — ► Eq. (6.12) is simply the expansion 

rate of a Universe filled with a perfect fluid p a . We stress that this solution is only valid 

in the neighborhood of the initial integration time to: later on the full numerical solution 

is mandatory. The mentioned features of the background solutions are illustrated in Fig. 

2, where the evolution of the common logarithms of the energy density are reported as a 

function of the common logarithm of the cosmic time coordinate 9 . As a side remark, we 

may note that in Fig. 2 (plot on the left-hand side) the first two lines at the top illustrate 

the behaviour of the function g a (t) (defined in Eq. (4.27)) and of its approximation (valid 

for large times) i.e. H/T. The function g & is always regular and never vanishes. This aspect 

is related to the occurrence that the contribution of the bulk viscosity becomes negligible, 

in practice, when the b-fluid dominates. 

9 In the remaining part of the article the logarithms will always be common logarithms, i.e. logarithms 
to base 10. 
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Figure 2: Evolution of the energy densities of the a- and b- fluid for different choices of the 
parameters in the case when the homogeneous bulk viscosity coefficient is parametrized as 
Ca = e y/Pa.- The thin lines always indicate the common logarithm of the energy density of 
the a-fluid, while the thick lines refer to the common logarithm of the energy density of the 
b-fluid (that dominates for t S> T ). 



It is now appropriate to analyse the evolution of C and Ca- Consider, to begin with, the 
case where e = 0. Following the discussion of [41] (see also [45]), the asymptotic value of ( is 
Cfinal — — <5r/6. This is exactly what is found numerically (see Fig. 3: lower curves in both 
plots at the left- and right-hand sides). In Fig. 3 the predicted asymptotic value is reported 
with a dashed line. As already remarked in [41] the accuracy of the estimate is controlled 
by the ratio T/Hq. This aspect is also illustrated in Fig. 3. In the plot on the left-hand 
side T/Hq = 10~ 3 , while on the right-hand side T/Hq = 0.1. Even if the asymptotic value 
of | C / <5r | is always well approximated by 1/6, the prediction is, comparatively, less accurate 
when T/Hq is larger. 

When e / and the discussion of [41] can be generalized. The authors of Ref. [41] notice 
that if is the solution of the system for <5p = 0, then the "forced" solution (corresponding 
to the case 5r ^ ) can be obtained as = — 5r- Assuming that, indeed, for large times 
both the curvature fluctuations £ and the Bardeen potential ^ go to a constant (as it is the 
case from the numerical analysis) , we have that the final asymptotic value of in the case 
where St ^ and e 7^ is given by 

xT/ _ [(3w a + 5) - 3V3e](w h + 1) - (3^ b + 5)[(w a + 1) - V3e] 

(3w h + 5)[{w a + 1) - V3e] 

Recall now that, when the a-fluid has already decayed, the relation between £ and (up 
to spatial gradients, i.e. for length-scales larger than the Hubble radius) is nothing but 

Cfinal = ~ o/ Wb li\ Cfinal- (6-15) 

3{w h + 1) 
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Figure 3: Evolution of curvature fluctuations in the case e = and for different values of 
the homogeneous component of the decay rate. 

From Eqs. (6.14) and (6.15), the final value of curvature fluctuations is then given by the 
following general equation: 

c f K + 1) (3w a + 5 - 3y/3e) - (3w h + 5) [(w a + l)-vge)] 1 

1 3(» b + i)[K + i)-V36] J r ' 1 " ) 

where by Cfinai we mean, as before, the final asymptotic value in the phase dominated by 
the decay products, i.e. the b-fluid. 

In the case w a = (leaving undetermined) Eq. (6.16) leads to 

[{Wb + 1)(5 - 3^/36) - {3w h + 5)(1 - y/3e) 



Co- 



■w b 



-S T . 



(6.17) 



3(u; b + 1)(1 - y^e) 

In Eq. (6.17) the subscript of £ means that the quantity computed at the right-hand side is 
the asymptotic value of £ in the case of the decay of w a = to ^ w a . The same notation 
will be used in the following equations. 

In the case = 1/3 (leaving u> a undetermined) Eq. (6.16) leads to 



Cw a ->i/3 



(1 - 3u> a ) + 3V3e 



+ 1- V3e 



Sr. 



(6.18) 



By setting w\, = 1/3 in Eq. (6.17) we will get the result valid for the case of a mixture 
of dust matter and radiation, i.e. 

1/1 + 3V3V 



If we set w a = 1 in Eq. (6.18), we deduce that 

1 2 - 3y/3e 



Sr- 



Cl->l/3 



6 2 - y/3e 



Sr- 



(6.19) 



(6.20) 
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Figure 4: The analytical results obtained from Eq. (6.16) are compared with the numerical 
calculation for the evolution of £ for the transition dust-radiation (left plot) and for the 
transition stiff-radiation. 



By setting e — > in Eq. (6.19) the usual result (valid in the absence of bulk viscosity) 
is reproduced. It is an amusing numerical coincidence that, in the limit e — > 0, Co->l/3 = 
-&-»i/3j i- e - IC0-4I/3I = IC1-+1/3I = 1/6. 

These results compare rather well with the numerical integration of the full system. In 
Fig. 4 we report the result for the numerical integration of the fluctuations in terms of 
the logarithm of | C/^r I ■ The plot at the left has to be compared with the prediction of Eq. 
(6.19) (appropriate for a transition dust-radiation transition), while the plot at the right has 
to be compared with Eq. (6.20) (appropriate for a stiff-radiation transition). In both plots 
the predictions of, respectively, Eqs. (6.19) and (6.20) is illustrated with a dashed line. For 
comparison, we also report the expectation obtainable in the limit e — > 0. The numerical 
correctness of Eq. (6.16) seems then numerically justified. In Fig. 5 the integration of the 
evolution equations of the fluctuations is reported for different values of e and 6r ■ 

In more general terms one can think of cases when ^ a = ep™. In the absence of a 
decay rate ( and for a single fluid) [7], if m > 1/2, there exist general deflationary solutions 
beginning in a quasi-de Sitter stage that evolve away from it for large times. Always in 
the absence of decay rate (and always in the single fluid case) [7], if m < 1/2 the solution 
starts expanding and reaches, for large times, a quasi-de Sitter solution. These statements 
hold, as stressed, for the case of a single fluid and in the absence of a decay rate. If there 
are many fluids, allowed to decay, the evolution depends upon the relative balance of the 
different parameters of the model. The main end of the present approach is to investigate 
the situations when the bulk viscosity coefficient does not dominate for large positive times. 
In this sense, solutions that do not lead to a radiation-dominated background for t > 1/T 
are less interesting for the present analysis (but may be relevant in related contexts). 
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Figure 5: The integration of the evolution equations of the fluctuations is reported for 
different mixtures of fluids. In the plot at the left 5p = 10 _6 ,at the right dp = 10~ 8 . 

In this perspective, even if m < 1/2, it could happen that radiation dominates very 
quickly if the decay rate is sufficiently large and e is sufficiently small. Consider, for instance, 
the case illustrated in Fig. 6, where we took £ a = epi /4 . If T/Hq = 0.1 and e = 10" 3 the 
evolution equations can be integrated. To be safe, the evolution equations of the fluctuations 
may be integrated directly in the off-diagonal gauge in terms of <p (not of Q). This choice 
has the advantage that the evolution of 4> is never singular even if H — > during a possible 
quasi-de Sitter stage. Also other cases of this rich system may be physically relevant. Those 
presented here were selected because they illustrate, in simple terms, the interplay between 
bulk viscous stresses and fluctuating decay rates. 

7 Concluding remarks 

Bulk viscous effects of interacting relativistic fluids have been investigated. These two effects 
were never considered simultaneously, to the best of our knowledge. The energy-momentum 
tensors of each of the fluids of the relativistic plasma have been parametrized as the sum 
of an inviscid contribution and of a bulk viscous correction. Different results have been 
obtained: 

• a fully gauge-invariant description of the evolution of the fluctuations of the geometry 
has been discussed in the presence of both a fluctuating decay rate and a fluctuating 
bulk viscosity coefficient; 

• gauge-dependent descriptions have been derived and exploited with particular atten- 
tion to the off-diagonal (or uniform-curvature) gauge. 
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Figure 6: The integration of the system is illustrated in the case £ a = ep a and £b = 0. 
The value 5r = 10 -8 has been used. In the plot on the left-hand side the evolution of the 
background energy density is reported on a double logarithmic scale. On the right-hand- 
side plot the evolution of C and £ a is reported for the same background values used in the 
left hand one. 

• if the bulk viscosity coefficient and the decay rate are allowed to have spatial variation 
on super-horizon scales, curvature fluctuations may be generated even if they are 
initially vanishing; 

• examples of the latter phenomenon were provided both analytically and numerically. 

As far as applications are concerned, two different situations may arise. In the first case 
the bulk viscosity is absent from the background, but its inhomogeneities may affect the 
evolution of the fluctuations with wavelengths larger than the Hubble radius. In the second 
case, the bulk viscosity does contribute to the background. It has been shown that situations 
can be envisaged where the contribution of the bulk viscosity coefficient becomes subleading 
after the decay of the viscous fluid of the mixture. In this case, the curvature fluctuations for 
length scales larger that the Hubble radius are determined by the interplay of the fluctuating 
viscosity and of the inhomogeneous decay rate. The final asymptotic value of the curvature 
fluctuations has been computed in specific numerical examples, where, for instance, the 
decay products are in the form of a radiation fluid. 
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Appendix: vector modes 



For completeness, some results concerning the evolution of vector modes in a mixture of 
viscous interacting fluids are reported in the present appendix. They are not central to the 
present analysis, but are logically connected to the considerations presented so far. In fact, 
one may wonder if effects similar to the ones discussed for the scalar modes may also arise 
in the case of vector modes. Consider then the vector fluctuations of the four-dimensional 
metric (3.4). They are parametrized by two (solenoidal) three-dimensional vectors Qi and 
Wi which have, overall, three independent components: 

c 2/-, r Oi Q % 

o v goi = -a Qi, d v g = 5-, 

Sygij = a 2 (diW 3 + djWi), 5 v g ij = -^(VW* + &W l ), (A.l) 

where 5 V denotes the vector fluctuation of the various entries of the perturbed metric. It is 
possible to select a gauge where S v gij is vanishing. This choice amounts to fixing two degrees 
of freedom of the perturbed metric. The vector fluctuations of the Einstein equations as 
well as the vector fluctuations of the equations describing the energy-momentum exchange 



between different fluids can be written as 

SvGl = 8ttG6 v T», (A.2) 

V^rr = -T5 v g ua u a ( Pa + p a ) - Tg ua S v u a ( Pa + p a ), (A.3) 

V M <5 v T b ^ = T5 v g ua u a (p a + p a ) + Tg va S v u a (p a + p a ) , (A.4) 



with the same notation as used in Eqs. (1.1), (1.2) and (2.4), but with the difference that, 
here, vector fluctuations are considered. 

Using the perturbed metric (A.l) and recalling that uq8 v u z = V 1 , the (Oi) and (ij) 
components of Eq. (A.2) imply: 

V 2 Qi = -WttG( P + V)a 2 Vi, (A.5) 

Qi + 2HQi = 0. (A.6) 

Equations (A.3) and (A.4) lead, instead, to the following pair of equations: 

[Vf } (p a + Pa)]' + m[vt ] (Pa + V a )} = -af(p a + Pa )V h (A.7) 
[Vf } (pb + Pb)]' + 4H[vf \p h + V h )\ = af(p a + Pa )Vi. (A.8) 

Considerations similar to the ones reported in Eqs. (3.42), (3.43) and (3.45) lead to the 
evolution equation of the total velocity field V«, i.e. 

[Vi(p a + Pa)]' + 4H[V,(p a + V a )} = 0, 

(p + V) V, = (pa + V & ) V, (a) + (p b + Vb) Vf } . (A.9) 
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From Eqs. (A.5)-(A.8) it can be easily argued that the spatial fluctuations of the 
bulk viscosity coefficient do not contribute, as expected, to the evolution of the vector 
fluctuations of the geometry; the situation is then very similar to the one discussed in the 
inviscid case in a number of studies that appeared in the past [55, 56] as well as more 
recently [57, 58, 59, 60]. The difference is, of course, that in the inviscid case the evolution 
equations of the background are different from those where a homogeneous bulk viscosity 
coeffcient is included. Finally, always for completeness, we should mention that the tensor 
modes, in some classes of bulk viscous solution, have been analysed in [61]. 
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